#pragma once
#ifndef MESH_HPP
#define MESH_HPP
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <thread>
#include <functional>
#include <mutex>
#include <ctime>
#include <cmath>
#include <iomanip>
#include <cfloat>
#include <climits>
#include "thread_manager.hpp"
#include "FileMapping.hpp"

#include <unordered_map>

using namespace std;
class CPgrid;
struct cpoint_hash;

#define HEXBOTTOM 0
#define HEXFRONT 1
#define HEXLEFT 2
#define HEXBACK 3
#define HEXRIGHT 4
#define HEXUP 5

//定义PI
#define M_PI_4_FAULTS  3.14159265358979323846
//定义epsilon
const double EPS = 1e-7;

//判断断层阈值角
#define FAULT_MIN_ANGLE 158
//断层阈值角的cos值
#define FAULT_MIN_COS cos(FAULT_MIN_ANGLE/180.0* M_PI_4_FAULTS)

//获取a向量与b向量夹角余弦值
#define get_cos_abs(a,b)  (a.x * b.x + a.y * b.y + a.z * b.z) / (sqrt(a.x * a.x + a.y * a.y + a.z * a.z) * sqrt(b.x * b.x + b.y * b.y + b.z * b.z))

#define get_cos_2d_abs(a,b)  (a[0] * b[0] + a[1] * b[1] ) / (sqrt(a[0] * a[0] + a[1] * a[1]) * sqrt(b[0] * b[0] + b[1] * b[1]))

#define cpoint_distance(a,b) ((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z))
// 通过z以及柱线，计算x y
// Qk=Qi+(Qj-Qi)*(Zk-Zi)/(Zj-Zi);
// Q=X/Y
#define getXY(Qi, Qj, Zi, Zj, Zk) (Qi + ((Qj - Qi) * (Zk - Zi) / (Zj - Zi)))

// 若x==a则x=b;
// 若x==b则x=a;
// 否则不变
#define local_change(x, a, b) (x = (x == a ? b : (x == b ? a : x)))

//判断x是否可能是数字
#define might_be_num(x) ((x>='0'&&x<='9')||(x=='.')||(x=='-'))

// 点
class cpoint
{
public:
    //点的坐标
    double x, y, z;
    //标记是否是断层点
    int is_fault = 0;
    //默认构造函数
    cpoint()
    {
        x = 0.0;
        y = 0.0;
        z = 0.0;
    }
    //普通构造函数
    cpoint(double a, double b, double c) : x(a), y(b), z(c) { ; }
    //复制构造函数
    cpoint(const cpoint& t)
    {
        x = t.x;
        y = t.y;
        z = t.z;
        is_fault = t.is_fault;
    }
    //重载 == 运算符
    bool operator==(const cpoint &t) const
    {
        return fabs(x-t.x)<EPS&&fabs(y - t.y)<EPS&&fabs(z - t.z)<EPS;
    }
    //重载 = 运算符
    cpoint operator=(cpoint t)
    {
        x = t.x;
        y = t.y;
        z = t.z;
        is_fault = t.is_fault;
        return *this;
    }
    //重载 - 运算符
    cpoint operator-(cpoint t)
    {
        cpoint res;
        res.x = x - t.x;
        res.y = y - t.y;
        res.z = z - t.z;
        return res;
    }

    //重载 + 运算符
    cpoint operator+(cpoint t)
    {
        cpoint res;
        res.x = x + t.x;
        res.y = y + t.y;
        res.z = z + t.z;
        return res;
    }

    //重载 / 运算符
    cpoint operator/(double t)
    {
        cpoint res;
        res.x = x/t;
        res.y = y/t;
        res.z = z/t;
        return res;
    }
    //重载 * 运算符
    cpoint operator*(double t)
    {
        cpoint res;
        res.x = x * t;
        res.y = y * t;
        res.z = z * t;
        return res;
    }
    cpoint product(cpoint& t) {
        cpoint k;
        k.x = y * t.z - t.y * z;
        k.y = t.x * z - t.z * x;
        k.z = x * t.y - t.x * y;
        return k;
    }

    bool on_the_segment(cpoint pnt1, cpoint pnt2) {
        double a, b, c;
        a = sqrt(cpoint_distance((*this), pnt1));
        b = sqrt(cpoint_distance((*this), pnt2));
        c = sqrt(cpoint_distance(pnt1, pnt2));
    }
};
struct cpoint_hash {
    size_t operator()(const cpoint& t) const {
        return hash<double>()(t.x) ^ hash<double>()(t.y) ^ hash<double>()(t.z);// hash<int>()(p.age);
    }
};
class cpline {
public:
    cpoint mypoint;
    cpoint dir;
    cpline() {
        mypoint = cpoint();
        dir = cpoint();
    }
    cpline(cpoint a, cpoint b, int mode = 0) {
        if (mode == 0) {
            mypoint = a;
            dir = b - a;
        }
        else {
            mypoint = a;
            dir = b;
        }
    }

    cpoint intersection(cpline a) {
        cpoint cof[2];
        cof[0] = a.mypoint;
        cof[1] = a.dir;
        double t1 = (cof[0].y - mypoint.y - dir.y * (cof[0].x - mypoint.x) / dir.x) / (dir.y * cof[1].x / dir.x - cof[1].y);
        if (fabs((dir.y * cof[1].x / dir.x - cof[1].y)) < EPS)
            cout << "wtf???\n";
        return cpoint(a.mypoint+a.dir*t1);
    }
};
// 点集
class pset
{
public:
    //点集数目
    int psize = 0;
    //点指针
    cpoint* pdata;
    pset() { pdata = nullptr; }
};
// 六面体单元
//该单元不存储自身id值，id值在CPgrid中体现：
//某一单元的id值等同于其在CPgrid::hexdata中的物理索引值
class cphex
{
public:

    cphex(){;}
    ~cphex() {
        //if(faults!=nullptr)
        //    delete[]faults;
    }
    cphex operator=(cphex t)
    {
        for (int i = 0; i < 8; i++)
            mpoint[i] = t.mpoint[i];
        has_faults = t.has_faults;
        for (int i = 0; i < has_faults; i++) {
            if (faults == nullptr)
                faults = new int[has_faults];
            faults[i] = t.faults[i];
        }
        child = t.child;
        center = t.center;
        stat = t.stat;
        ext = t.ext;
        valid = t.valid;
        poissonr = t.poissonr;
        youngmode = t.youngmode;
        max_edge_len = t.max_edge_len;
        rockden=t.rockden;
        poro=t.poro;
        return *this;
    }

    // 仅仅存储点的编号
    // 存储顺序以及空间相对位置如下
    // bottom:   3---2    top:   7---6
    //           |   |           |   |
    //           |   |           |   |
    //           0---1           4---5
    int mpoint[8] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN, INT_MIN, INT_MIN, INT_MIN, INT_MIN};

    double max_edge_len = DBL_MIN;
    //记录该单元点中属于断层处的点数目
    int has_faults = 0;
    //记录断层点的id
    int* faults=nullptr;

    cpoint center;
    //判断该单元是否为断层附近单元，同时在faults中记录本单元中的断层点
    void cal_fault(cpoint* pnts, int mode = 0, int *m_faults =nullptr);
    void cal_center(cpoint* pnts) {
        center.x = 0.0;
        center.y = 0.0;
        center.z = 0.0;
        for (int i = 0; i < 8; i++) {
            center.x += pnts[mpoint[i]].x;
            center.y += pnts[mpoint[i]].y;
            center.z += pnts[mpoint[i]].z;
        }
        center.x /= 8.0;
        center.y /= 8.0;
        center.z /= 8.0;
    };
    void cal_size(cpoint* pnts) {
        cpoint temp;
        max_edge_len = DBL_MIN;
        for (int i = 0; i < 12; i++) {
            switch (i) {
            case 0:
                temp = pnts[mpoint[0]] - pnts[mpoint[1]];
                break;
            case 1:
                temp = pnts[mpoint[0]] - pnts[mpoint[3]];
                break;
            case 2:
                temp = pnts[mpoint[2]] - pnts[mpoint[3]];
                break;
            case 3:
                temp = pnts[mpoint[1]] - pnts[mpoint[2]];
                break;
            case 4:
                temp = pnts[mpoint[0]] - pnts[mpoint[4]];
                break;
            case 5:
                temp = pnts[mpoint[3]] - pnts[mpoint[7]];
                break;
            case 6:
                temp = pnts[mpoint[2]] - pnts[mpoint[6]];
                break;
            case 7:
                temp = pnts[mpoint[1]] - pnts[mpoint[5]];
                break;
            case 8:
                temp = pnts[mpoint[4]] - pnts[mpoint[5]];
                break;
            case 9:
                temp = pnts[mpoint[4]] - pnts[mpoint[7]];
                break;
            case 10:
                temp = pnts[mpoint[6]] - pnts[mpoint[7]];
                break;
            case 11:
                temp = pnts[mpoint[5]] - pnts[mpoint[6]];
                break;
            default:
                break;
            }
            max_edge_len = max(max_edge_len,sqrt(temp.x * temp.x + temp.y * temp.y + temp.z * temp.z));
        }
    }

    int stat = -1;
    cphex **neighbors=nullptr;
    int child_face[6];
    int** child_tetra = nullptr;
    double** tetra_property = nullptr;
    cphex* child = nullptr;
    int n_child = 0;
    void optimize(vector<cpoint> &pdata, unordered_map<cpoint, int, cpoint_hash> &point_map) {
        int m_faults[8] = { 0 };
        cal_fault(&pdata[0],1,m_faults);
        int badid = -1;

        for (int i = 0; i < 8; i++) {
            if (m_faults[i] !=0) {
                badid = i;
                break;
            }
        }

        int half = 0;
        int up = 0;
        int down = 0;
        for (int i = 0; i < 4; i++) {
            up += m_faults[i];
            down += m_faults[i + 4];
        }

        if ((up == 0 && down != 0) || (up != 0 && down == 0)) {
            if (up == 0 && down != 0)
                half = 1;
            else
                half = 2;
        }

        if (badid != -1) {
            cpoint pnts[8];
            for (int i = 0; i < 4; i++) {
                pnts[i] = pdata[mpoint[(badid + i) % 4]];
                pnts[i + 4] = pdata[mpoint[(badid + i) % 4 + 4]];
            }
            cpoint mid[4];
            cpoint new_hex[8];
            cpoint line01, line02, line03;
            double angle102, angle203, angle103;
            int new_id[8];
            int topo_id[8];
            for (int i = 0; i < 4; i++) {
                topo_id[i] = mpoint[(badid + i) % 4];
            }

            for (int i = 4; i < 8; i++) {
                topo_id[i] = mpoint[(badid + i) % 4 + 4];
            }

            mid[0] = (pnts[1] + pnts[2]) / 2;
            mid[1] = (pnts[3] + pnts[2]) / 2;
            mid[2] = (pnts[6] + pnts[5]) / 2;
            mid[3] = (pnts[6] + pnts[7]) / 2;
            for (int i = 0; i < 4; i++) {
                double ang0, ang1,ang2;
                cpoint a, b;
                cpoint p0, p1, p2;
                switch (i)
                {
                case 0: {
                    a = pnts[1];
                    b = pnts[2];
                    p0 = pnts[2];
                    p1 = pnts[3];
                    p2 = pnts[0];
                    break;
                }
                case 1: {
                    a = pnts[3];
                    b = pnts[2];
                    p0 = pnts[2];
                    p1 = pnts[1];
                    p2 = pnts[0];
                    break;
                }
                case 2: {
                    a = pnts[5];
                    b = pnts[6];
                    p0 = pnts[6];
                    p1 = pnts[7];
                    p2 = pnts[4];
                    break;
                }
                case 3: {
                    a = pnts[7];
                    b = pnts[6];
                    p0 = pnts[6];
                    p1 = pnts[5];
                    p2 = pnts[4];
                    break;
                }
                default:
                    break;
                }
                mid[i] = (a + b) / 2;
                ang0 = acos(get_cos_abs((p0 - p1), (mid[i] - p1)));
                ang1 = acos(get_cos_abs((p2 - p1), (mid[i] - p1)));
                ang2 = acos(get_cos_abs((p2 - p1), (p0 - p1)));

                while ((fabs(ang0 - ang1) > 0.05 || (ang0 + ang1) > ang2+0.05)&&!(mid[i]==a||mid[i]==b)) {
                    if (ang0>ang1) {
                        a = mid[i];
                        mid[i] = (a + b) / 2;
                    }
                    else {
                        b = mid[i];
                        mid[i] = (a + b) / 2;
                    }
                    ang0 = acos(get_cos_abs((p0 - p1), (mid[i] - p1)));
                    ang1 = acos(get_cos_abs((p2 - p1), (mid[i] - p1)));
                }
            }

            double ratio_to_bad = 0.9;
            double ratio_to_corner = 0.3;
            double ratio_left = 0.8;
            double ratio_right = 0.8;
            new_hex[0] = (pnts[0] - pnts[2]) * ratio_to_bad + pnts[2];
            new_hex[1] = (pnts[1] - mid[1]) * ratio_right + mid[1];
            new_hex[2] = (pnts[0] - pnts[2]) * ratio_to_corner + pnts[2];
            new_hex[3] = (pnts[3] - mid[0]) * ratio_left + mid[0];

            line01 = new_hex[1] - new_hex[0];
            line02 = new_hex[2] - new_hex[0];
            line03 = new_hex[3] - new_hex[0];
            angle103 = get_cos_abs(line01, line03);
            angle102 = get_cos_abs(line01, line02);
            angle203 = get_cos_abs(line02, line03);

            while ((acos(get_cos_abs(line02, (pnts[1] - new_hex[0]))) >= (M_PI_4_FAULTS / 4) - 0.05) 
                &&(acos(angle102) >= (M_PI_4_FAULTS / 4) - 0.05)) {
                
                ratio_right -= 0.001;
                new_hex[1] = (pnts[1] - mid[1]) * ratio_right + mid[1];
                line01 = new_hex[1] - new_hex[0];
                angle102 = get_cos_abs(line01, line02);
            }

            while ((acos(get_cos_abs(line02,(pnts[3]-new_hex[0]))) >= (M_PI_4_FAULTS / 4) - 0.05)
                && (acos(angle203) >= (M_PI_4_FAULTS / 4) - 0.05)) {
                ratio_left -= 0.001;
                new_hex[3] = (pnts[3] - mid[0]) * ratio_left + mid[0];
                line03 = new_hex[3] - new_hex[0];
                angle203 = get_cos_abs(line02, line03);
            }

            while (angle103 < FAULT_MIN_COS || acos(angle203) + acos(angle102) > M_PI_4_FAULTS) {
                if (angle102 > angle203) {
                    if ((acos(angle203) > M_PI_4_FAULTS / 3) && cpoint_distance(new_hex[0], new_hex[1]) < cpoint_distance(new_hex[0], new_hex[3])) {
                        ratio_left -= 0.01;
                        new_hex[3] = (pnts[3] - mid[0]) * ratio_left + mid[0];
                        line03 = new_hex[3] - new_hex[0];
                    }
                    else {
                        ratio_right -= 0.01;
                        new_hex[1] = (pnts[1] - mid[1]) * ratio_right + mid[1];
                        line01 = new_hex[1] - new_hex[0];
                    }
                }
                else {
                    if ((acos(angle102) > M_PI_4_FAULTS / 3) && cpoint_distance(new_hex[0], new_hex[1]) > cpoint_distance(new_hex[0], new_hex[3])) {
                        ratio_right -= 0.01;
                        new_hex[1] = (pnts[1] - mid[1]) * ratio_right + mid[1];
                        line01 = new_hex[1] - new_hex[0];
                    }
                    else {
                        ratio_left -= 0.01;
                        new_hex[3] = (pnts[3] - mid[0]) * ratio_left + mid[0];
                        line03 = new_hex[3] - new_hex[0];
                    }
                }
                angle102 = get_cos_abs(line01, line02);
                angle203 = get_cos_abs(line02, line03);
                angle103 = get_cos_abs(line01, line03);
            }

            ratio_to_bad = 0.9;
            ratio_to_corner = 0.3;
            ratio_left = 0.8;
            ratio_right = 0.8;

            new_hex[4] = (pnts[4] - pnts[6]) * ratio_to_bad + pnts[6];
            new_hex[5] = (pnts[5] - mid[3]) * ratio_right + mid[3];
            new_hex[6] = (pnts[4] - pnts[6]) * ratio_to_corner + pnts[6];
            new_hex[7] = (pnts[7] - mid[2]) * ratio_left + mid[2];

            line01 = new_hex[5] - new_hex[4];
            line02 = new_hex[6] - new_hex[4];
            line03 = new_hex[7] - new_hex[4];
            angle103 = get_cos_abs(line01, line03);
            angle102 = get_cos_abs(line01, line02);
            angle203 = get_cos_abs(line02, line03);
            
            while ((acos(get_cos_abs(line02, (pnts[5] - new_hex[4]))) >= (M_PI_4_FAULTS / 4) - 0.05)
                && (acos(angle102) >= (M_PI_4_FAULTS / 4) - 0.05)) {

                ratio_right -= 0.001;
                new_hex[5] = (pnts[5] - mid[3]) * ratio_right + mid[3];
                line01 = new_hex[5] - new_hex[4];
                angle102 = get_cos_abs(line01, line02);
            }

            while ((acos(get_cos_abs(line02, (pnts[7] - new_hex[4]))) >= (M_PI_4_FAULTS / 3) - 0.05)
                && (acos(angle203) >= (M_PI_4_FAULTS / 4) - 0.05)) {

                ratio_left -= 0.001;
                new_hex[7] = (pnts[7] - mid[2]) * ratio_left + mid[2];
                line03 = new_hex[7] - new_hex[4];
                angle203 = get_cos_abs(line02, line03);
            }
        

            while (angle103 < FAULT_MIN_COS || acos(angle203) + acos(angle102) > M_PI_4_FAULTS) {
                if (angle102 > angle203) {
                    if ((acos(angle203) > M_PI_4_FAULTS / 3) && cpoint_distance(new_hex[4], new_hex[5]) < cpoint_distance(new_hex[4], new_hex[7])) {
                        ratio_left -= 0.01;
                        new_hex[7] = (pnts[7] - mid[2]) * ratio_left + mid[2];
                        line03 = new_hex[7] - new_hex[4];
                    }
                    else {
                        ratio_right -= 0.01;
                        new_hex[5] = (pnts[5] - mid[3]) * ratio_right + mid[3];
                        line01 = new_hex[5] - new_hex[4];
                    }
                }
                else {
                    if ((acos(angle102) > M_PI_4_FAULTS / 3) && cpoint_distance(new_hex[4], new_hex[5]) > cpoint_distance(new_hex[4], new_hex[7])) {
                        ratio_right -= 0.01;
                        new_hex[5] = (pnts[5] - mid[3]) * ratio_right + mid[3];
                        line01 = new_hex[5] - new_hex[4];
                    }
                    else {
                        ratio_left -= 0.01;
                        new_hex[7] = (pnts[7] - mid[2]) * ratio_left + mid[2];
                        line03 = new_hex[7] - new_hex[4];
                    }
                }
                angle102 = get_cos_abs(line01, line02);
                angle203 = get_cos_abs(line02, line03);
                angle103 = get_cos_abs(line01, line03);
            }

            if (half != 0)
                n_child = 6;
            else
                n_child = 5;
            
            child = new cphex[n_child];

            unordered_map<cpoint, int, cpoint_hash>::iterator it;

            if (half != 0) {
                if (half == 1) {
                    for (int i = 0; i < 4; i++) {
                        new_hex[i] = (new_hex[i + 4] - new_hex[i]) * 0.3 + new_hex[i];
                    }
                }
                else if (half == 2) {
                    for (int i = 4; i < 8; i++) {
                        new_hex[i] = (new_hex[i-4] - new_hex[i]) * 0.3 + new_hex[i];
                    }
                }
            }

            for (int i = 0; i < 8; i++) {
                it = point_map.find(new_hex[i]);
                if (it == point_map.end()) {
                    child[4].mpoint[i] = pdata.size();
                    new_id[i] = pdata.size();
                    point_map.insert({ new_hex[i], pdata.size() });
                    pdata.push_back(new_hex[i]);
                }
                else {
                    child[4].mpoint[i] = it->second;
                    new_id[i] = it->second;
                }
            }

            if (half != 0) {
                if (half == 1) {
                    child[5].mpoint[0] = topo_id[0];
                    child[5].mpoint[1] = topo_id[1];
                    child[5].mpoint[2] = topo_id[2];
                    child[5].mpoint[3] = topo_id[3];
                    child[5].mpoint[4] = new_id[0];
                    child[5].mpoint[5] = new_id[1];
                    child[5].mpoint[6] = new_id[2];
                    child[5].mpoint[7] = new_id[3];
                }
                else if (half == 2) {
                    child[5].mpoint[0] = new_id[4];
                    child[5].mpoint[1] = new_id[5];
                    child[5].mpoint[2] = new_id[6];
                    child[5].mpoint[3] = new_id[7];
                    child[5].mpoint[4] = topo_id[4];
                    child[5].mpoint[5] = topo_id[5];
                    child[5].mpoint[6] = topo_id[6];
                    child[5].mpoint[7] = topo_id[7];
                }
            }

            for (int i = 0; i < 4; i++) {
                child[i].mpoint[0] = topo_id[i];
                child[i].mpoint[1] = topo_id[(i + 1) % 4];;
                child[i].mpoint[2] = new_id[(i + 1) % 4];
                child[i].mpoint[3] = new_id[i];
                child[i].mpoint[4] = topo_id[i + 4];
                child[i].mpoint[5] = topo_id[(i + 1) % 4 + 4];
                child[i].mpoint[6] = new_id[(i + 1) % 4 + 4];
                child[i].mpoint[7] = new_id[4+i];
            }
            for (int i = 0; i < n_child; i++) {
                child[i].ext = ext;
            }
        }
    }

    void optimize2tetra(vector<cpoint>& pdata,unordered_map<cpoint, int, cpoint_hash>& point_map) {
        int m_faults[8] = { 0 };
        cal_fault(&pdata[0], 1, m_faults);
        int badid = -1;

        int c = 0;
        for (int i = 0; i < 8; i++) {
            if (m_faults[i] != 0) {
                badid = i;
                ++c;
                break;
            }
        }
        
        if (badid != -1) {
            int topo_hexid[8] = { 0 };

            for (int q = 0; q < 4; q++) {
                topo_hexid[q] = mpoint[(q + badid) % 4];
                topo_hexid[q + 4] = mpoint[(q + badid) % 4 + 4];
            }
            child_tetra = new int* [6];
            tetra_property = new double* [6];

            for (int i = 0; i < 6; i++) {
                child_tetra[i] = new int[4];
                tetra_property[i] = new double[4];
                for (int j = 0; j < 4; j++) {
                    tetra_property[i][j] = 0.0;
                }
            }

            child_tetra[0][0] = topo_hexid[4];
            child_tetra[0][1] = topo_hexid[7];
            child_tetra[0][2] = topo_hexid[6];
            child_tetra[0][3] = topo_hexid[3];

            child_tetra[1][0] = topo_hexid[4];
            child_tetra[1][1] = topo_hexid[0];
            child_tetra[1][2] = topo_hexid[3];
            child_tetra[1][3] = topo_hexid[2];

            child_tetra[2][0] = topo_hexid[4];
            child_tetra[2][1] = topo_hexid[2];
            child_tetra[2][2] = topo_hexid[3];
            child_tetra[2][3] = topo_hexid[6];

            child_tetra[3][0] = topo_hexid[4];
            child_tetra[3][1] = topo_hexid[6];
            child_tetra[3][2] = topo_hexid[5];
            child_tetra[3][3] = topo_hexid[1];

            child_tetra[4][0] = topo_hexid[4];
            child_tetra[4][1] = topo_hexid[1];
            child_tetra[4][2] = topo_hexid[0];
            child_tetra[4][3] = topo_hexid[2];

            child_tetra[5][0] = topo_hexid[4];
            child_tetra[5][1] = topo_hexid[1];
            child_tetra[5][2] = topo_hexid[2];
            child_tetra[5][3] = topo_hexid[6];

        }

    }

    // 标记该单元是否有效
    int valid = 1;
    //标记是否为延拓单元
    int ext = 1;

    //泊松比
    double poissonr = 0.0;
    //杨氏模量
    double youngmode = 0.0;

    double rockden=0.0;

    double poro=0.0;
    // 0 bottom
    // 1 front
    // 2 left
    // 3 back
    // 4 right
    // 5 up
    int split_face_stat[6][4] = {
        {0,1,2,3},
        {0,1,4,5},
        {0,3,4,7},
        {3,2,7,6},
        {1,2,5,6},
        {4,5,6,7}};

    // bottom        side             top
    //    -- 2 --      5-------6      -- 10--
    //   |       |     |       |     |       |
    //   1       3     |       |     9       11
    //   |       |     |       |     |       |
    //    -- 0 --      4-------7      -- 8 --
    int split_edge_stat[12][2] = {
        {0,1},
        {0,3},
        {2,3},
        {1,2},
        {0,4},
        {3,7},
        {2,6},
        {1,5},
        {4,5},
        {4,7},
        {7,6},
        {6,5} };

    void split(cphex *&sub, vector<cpoint>& pdata,int &psize,unordered_map<cpoint,int,cpoint_hash> &point_map);

    // pos belike  1 root -> 8 subs
    // down                 UP
    //  ------ ------        ------ ------
    // |      |      |      |      |      |
    // |   1  |   2  |      |   5  |   6  |
    // |      |      |      |      |      |
    //  ------ ------        ------ ------
    // |      |      |      |      |      |
    // |   0  |   3  |      |   4  |   7  |
    // |      |      |      |      |      |
    //  ------ ------        ------ ------

    // pos belike  1 root -> 8 subs
    // down                 UP
    //  ------ ------        ------ ------
    // |      |      |      |      |      |
    // |   2  |   3  |      |   6  |   7  |
    // |      |      |      |      |      |
    //  ------ ------        ------ ------
    // |      |      |      |      |      |
    // |   0  |   1  |      |   4  |   5  |
    // |      |      |      |      |      |
    //  ------ ------        ------ ------
    void sub_init(int new_stat,int pos,int *root_mpoint,int *id_point_on_face,int *id_point_on_edge,int centerid);
};


// 角点网格
class CPgrid : public pset
{
public:
    cpoint* get_pdata() { return pdata; }
    void init_neighbors(int hexid);

    vector<double> mark_base(const char* filename, int n_thread = 0);

    void mark_fault(const char *file);
    bool check_fault(cpoint& pnt);

    void add_cpoint(cpoint* t, int size);

    double coordinate_offset[2];

    // coordinate_trig_val[0] for cos
    // coordinate_trig_val[1] for sin
    double coordinate_trig_val[2];

    ~CPgrid() {
        //if(pdata!=NULL)
        //    delete[]pdata;
        //if(hexdata)
        //    delete[]hexdata;
        //if(ijline)
        //    delete[]ijline;
    }
    void clear() {
        delete[]pdata;
        delete[]hexdata;
        delete[]ijline;
    }

    cpoint cpgrid_bbox[2] = { cpoint(DBL_MAX,DBL_MAX,DBL_MAX),cpoint(-1.0*DBL_MAX,-1.0*DBL_MAX,-1.0*DBL_MAX) };

    double max_edge_len = DBL_MIN;
    // 去重后点的数目
    int u_psize = 0;
    // 单元数目
    int hexsize = 0;
    // 有效单元数目
    int hexsize_valid = 0;

    //六面体单元指针
    cphex *hexdata = nullptr;


    //初始化cpgrid
    void init(const char* filename, int n_thread = 0);
    void init(const char* input_file, const char* base, int n_thread = 0);
    //读入属性文件
    void read_properties(char* filename);

    //通过单元的ijk各方向索引，获取其指针
    inline cphex* get_hex(int x, int y, int z);
    //通过id获取索引的
    void get_hex_ijk(int id,int &x,int &y,int &z);

    int get_hex_id(int x, int y, int z);
    //通过给定id来获取左右前后上下邻接单元id
    void get_neighbor(int id, int* neighbors);

    vector<int> get_neighbor_range(vector<int> &gid,int incr=0);
    //准备adt建树所需数据接口
    void pre_ele_bbox_info(double* ele_bbox);
    void pre_adt_info(int* ele_list, double* nodelist);

    //输出vtk文件
    void trans2vtk(const char* out);

    //输出赛题一要求属性文件
    void out_property(char* out);

//private:
//
//文件名
    char filename[1024];

    //线程管理
    thread_manager threads;
    //子线程数
    int threads_num = 0;
    //总线程数=子线程数+1
    int all_threads = 0;

    //ijk方向尺寸
    int i = 0;
    int j = 0;
    int k = 0;

    //柱线的数量
    int ijsize = 0;
    //存储柱线两端点指针，
    cpoint* ijline = nullptr;
    vector<vector<cpoint>> fault_ij;
    double fault_bbox[2][3];
    //该变量用来标志是否读入文件是大文件；后续考虑优化大文件io
    bool bigfile = 0;
    //用于存储文件中zcorn字段所有double值
    double* ori_data=nullptr;

    FileMapping file_mmap;
    //各线程读入有效数据的指针
    double** threads_data = nullptr;
    //各线程读入有效数据统计值
    int* threads_data_valid = nullptr;
    //threads_data的size
    int* threads_data_size = nullptr;
    //zcorn字段在文件中出现位置
    streampos z_start_infile = 0;
    //zcorn字段在文件中结束位置
    streampos z_end_infile = 0;
    //小文件相关
    //标记zcorn字段在char*中的开始位置
    int z_start_inchar = 0;
    //标记zcorn字段在char*中的结束位置
    int z_end_inchar = 0;
    //用于存储文件内容
    char* file_content=nullptr;
    //文件长度
    long long file_len = 0;

    // 去重，构建的拓扑与要求的网格顶底面顺序一致
    // 0为主线程，其他为子线程
    void deal_zcorn(char* filename, int tid = 0);

    //检查给定线程数
    void thread_num_check();

    //读入zcorn
    void read_zcorn(int tid = 0);

    //处理除开zcorn外信息
    void cpgrid_info(char* filename);

    //判断hexid中的No.pnt_stat 点，是否与相邻其他单元中的可能点重合
    int point_exist(int hexid, int pnt_stat, double z, bool& mismatch);

    //单元中点的状态
    int s_stat_len[8] = {7, 6, 5, 4, 7, 6, 5, 4};
    //不同状态对应的，相关单元其ijk方向索引的偏移量，以及该点在相关单元中的位置
    int s_stat[8][7][4] = {
        {{-1, -1, -1, 7}, {0, -1, -1, 6}, {-1, 0, -1, 5}, {0, 0, -1, 4}, {-1, -1, 0, 3}, {0, -1, 0, 2}, {-1, 0, 0, 1}},
        {{0, -1, -1, 7}, {1, -1, -1, 6}, {0, 0, -1, 5}, {1, 0, -1, 4}, {0, -1, 0, 3}, {1, -1, 0, 2}},
        {{-1, 0, -1, 7}, {0, 0, -1, 6}, {-1, 1, -1, 5}, {0, 1, -1, 4}, {-1, 0, 0, 3}},
        {{0, 0, -1, 7}, {1, 0, -1, 6}, {0, 1, -1, 5}, {1, 1, -1, 4}},
        {{-1, -1, 0, 7}, {0, -1, 0, 6}, {-1, 0, 0, 5}, {-1, -1, 1, 3}, {0, -1, 1, 2}, {-1, 0, 1, 1}, {0, 0, 1, 0}},
        {{0, -1, 0, 7}, {1, -1, 0, 6}, {0, -1, 1, 3}, {1, -1, 1, 2}, {0, 0, 1, 1}, {1, 0, 1, 0}},
        {{-1, 0, 0, 7}, {-1, 0, 1, 3}, {0, 0, 1, 2}, {-1, 1, 1, 1}, {0, 1, 1, 0}},
        {{0, 0, 1, 3}, {1, 0, 1, 2}, {0, 1, 1, 1}, {1, 1, 1, 0}}};

    //根据传入的num类型，从file_content中获取一个数字，同时变更表明位置的ptr
    //仅支持int和double
    template<class type>
    inline type get_a_num(const char*file_content,long long &ptr,type num) {
        type res = num;
        //跳过非数字相关字符
        while (file_content[ptr]!='\0'&&!might_be_num(file_content[ptr]))
        {
            ++ptr;
        }
        //获取数字
        if (typeid(num)==typeid(int)) {
            res = atoi(file_content + ptr);
        }
        else if (typeid(num) == typeid(double)) {
            res = atof(file_content + ptr);
        }
        //跳过读入的数字
        while (might_be_num(file_content[ptr]))
        {
            ++ptr;
        }
        //返回
        return res;
    }


};


//拓扑六面体
class topo_hex {
public:
    //存储各方向的实际点
    int** direct[6];
    int *up[4];
    int *front[4];
    int *left[4];
    int *back[4];
    int *right[4];
    int *bottom[4];
    int *mpoint[8];
    //实际六面体
    cphex* host;
    topo_hex() { ; }

    topo_hex(cphex* t) {
        init(t);
    }

    //根据实际六面体初始化拓扑六面体
    void init(cphex *t) {
        host = t;
        direct[HEXUP] = up;
        direct[HEXBACK] = back;
        direct[HEXBOTTOM] = bottom;
        direct[HEXFRONT] = front;
        direct[HEXLEFT] = left;
        direct[HEXRIGHT] = right;

        up[0] = &t->mpoint[4];
        up[1] = &t->mpoint[5];
        up[2] = &t->mpoint[6];
        up[3] = &t->mpoint[7];
        bottom[0] = &t->mpoint[0];
        bottom[1] = &t->mpoint[1];
        bottom[2] = &t->mpoint[2];
        bottom[3] = &t->mpoint[3];

        front[0] = &t->mpoint[0];
        front[1] = &t->mpoint[1];
        front[2] = &t->mpoint[4];
        front[3] = &t->mpoint[5];
        back[0] = &t->mpoint[3];
        back[1] = &t->mpoint[2];
        back[2] = &t->mpoint[7];
        back[3] = &t->mpoint[6];

        left[0] = &t->mpoint[0];
        left[1] = &t->mpoint[3];
        left[2] = &t->mpoint[4];
        left[3] = &t->mpoint[7];
        right[0] = &t->mpoint[1];
        right[1] = &t->mpoint[2];
        right[2] = &t->mpoint[5];
        right[3] = &t->mpoint[6];

        mpoint[0] = bottom[0];
        mpoint[1] = bottom[1];
        mpoint[2] = bottom[2];
        mpoint[3] = bottom[3];

        mpoint[4] = up[0];
        mpoint[5] = up[1];
        mpoint[6] = up[2];
        mpoint[7] = up[3];
    }

    //定义当前拓扑六面体的邻居
    void is_neighbor_of(topo_hex& who,int direct) {
        int** t;
        int** neighbor_face;
        //根据相邻方向，赋值邻接面
        switch (direct) {
        case 0: {
            t = up;
            neighbor_face = who.bottom;
            break;
        }
        case 1: {
            t = back;
            neighbor_face = who.front;
            break;
        }
        case 2: {
            t = right;
            neighbor_face = who.left;
            break;
        }
        case 3: {
            t = front;
            neighbor_face = who.back;
            break;
        }
        case 4: {
            t = left;
            neighbor_face = who.right;
            break;
        }
        case 5: {
            t = bottom;
            neighbor_face = who.up;
            break;
        }
        }
        //建立邻接拓扑关系
        for (int i = 0; i < 4; i++) {
            *t[i] = *neighbor_face[i];
        }
    }
};

//有限元网格
class finiteGrid{
public:
    //点集
    vector<cpoint> pdata;
    //实际六面体
    vector<cphex>  hexdata;
    cpoint grid_bbox[2];
    int hexsize;
    int psize;

    //0 up ;; 1 front ;;2 left ;;3 back ;;4 right ;;5 down
    //点集哈希表，用于去重
    unordered_map<cpoint,int,cpoint_hash> point_map;
    //单元最长边
    double max_edge_len =DBL_MIN;
    
    //加密次数
    int ref_times=0;
    //未加密前，原始网格各方向层数
    int ori_ijk[3];
    //储层ijk方向层数
    int topo_i;
    int topo_j;
    int topo_k;
    //整体网格ijk方向层数
    int i_withext = 0;
    int j_withext = 0;
    int k_withext = 0;
    //实际延拓距离
    double leni_ext = 0.0;
    double lenj_ext = 0.0;
    double lenk_up_ext = 0.0;
    double lenk_bot_ext = 0.0;
    //含延拓网格中储层ijk方向索引起始值
    int reservior_relative_i_start = 0;
    int reservior_relative_j_start = 0;
    int reservior_relative_k_start = 0;
    //含延拓网格中储层ijk方向索引结束值
    int reservior_relative_i_end = 0;
    int reservior_relative_j_end = 0;
    int reservior_relative_k_end = 0;

    //延拓层生长比率
    //0 ->i;
    //1 ->j;
    //2 ->kup;
    //3 ->kdown;
    double height_ratio[4] = { 0.0,0.0,0.0,0.0 };
    //ijk方向延拓第一层高度
    double i_ori_height = 0.0;
    double j_ori_height = 0.0;
    double k_up_ori_height = 0.0;
    double k_bot_ori_height = 0.0;
    finiteGrid(){;}
    finiteGrid(CPgrid &ori_grid,int required_size){
        init(ori_grid,required_size);
    }

    //使用CPgrid初始化finiteGrid，同时设定需求尺寸
    void init(CPgrid &ori_grid,double required_size){
        //使用CPgrid相关数据初始化finiteGrid
        psize  =ori_grid.u_psize;
        hexsize=ori_grid.hexsize;
        hexdata.resize(hexsize);
        pdata.resize(psize);
        max_edge_len = ori_grid.max_edge_len;
        grid_bbox[0] = ori_grid.cpgrid_bbox[0];
        grid_bbox[1] = ori_grid.cpgrid_bbox[1];
        memcpy(&hexdata[0], ori_grid.hexdata, sizeof(cphex) * ori_grid.hexsize);
        memcpy(&pdata[0], ori_grid.pdata, sizeof(cpoint) * ori_grid.u_psize);
        //初始化哈希表
        for (int i = 0; i < psize; i++) {
            point_map.insert({pdata[i],i});
        }

        topo_i=ori_grid.i;
        topo_j=ori_grid.j;
        topo_k=ori_grid.k;
        
        ori_ijk[0]=ori_grid.i;
        ori_ijk[1]=ori_grid.j;
        ori_ijk[2]=ori_grid.k;

        reservior_relative_i_end = reservior_relative_i_start + topo_i;
        reservior_relative_j_end = reservior_relative_j_start + topo_j;
        reservior_relative_k_end = reservior_relative_k_start + topo_k;

        //加密操作
        refine(required_size);
    }

    //根据required_size和max_edge_len大小加密网格
    void refine(double required_size){
        //临时变量，用于存储六面体
        vector<cphex> t_hexdata;
        //当前ijk方向尺寸
        int cur_i = topo_i;
        int cur_j = topo_j;
        int cur_k = topo_k;
        int subid;
        int sub_ijk[3];
        int root_ijk[3];
        cphex* subhex = new cphex[8];

        double next_max_edge_len=DBL_MIN;
        while(max_edge_len >required_size){
            next_max_edge_len = DBL_MIN;
            t_hexdata.resize(cur_i * cur_j * cur_k * 8);
            //对所有六面体在三方向上二分加密
            for (int i = 0; i < hexdata.size(); i++) {
                hexdata[i].stat = i;
                hexdata[i].split(subhex,pdata,psize,point_map);
                //根据当前ijk尺寸，以及六面体线性索引i，获取其当前的ijk坐标
                get_ijk(i, cur_i, cur_j, cur_k, root_ijk[0], root_ijk[1], root_ijk[2]);
                for (int pos = 0; pos < 8; pos++) {
                    //获取加密后的子六面体ijk坐标
                    sub_ijk[0] = root_ijk[0] * 2 +  pos % 2;
                    sub_ijk[1] = root_ijk[1] * 2 + (pos % 4) / 2;
                    sub_ijk[2] = root_ijk[2] * 2 + pos / 4;
                    //根据子六面体ijk坐标，获取其线性索引，并存储至对应位置
                    subid = get_hex_id(sub_ijk[0], sub_ijk[1], sub_ijk[2], cur_i * 2, cur_j * 2, cur_k * 2);

                    subhex[pos].cal_size(&(pdata[0]));
                    t_hexdata[subid] = subhex[pos];
                    //更新最长边边长
                    next_max_edge_len = max(next_max_edge_len, subhex[pos].max_edge_len);
                }
            }

            max_edge_len = next_max_edge_len;
            //将新生成的六面体数据存储至hexdata中
            hexdata.clear();
            hexdata.assign(t_hexdata.begin(),t_hexdata.end());
            //更新储层ijk方向尺寸
            cur_i *= 2;
            cur_j *= 2;
            cur_k *= 2;
            ++ref_times;

            int m_faults[8] = { 0 };
            int ti, tj, tk;
            //对低质量单元进行保持形状操作
            //低质量单元直接进行二分加密，次数过多后，会产生边线切割顶点的情况
            //下面代码段对低质量单元进行了保持形状操作，避免出现如上情况
            {
                int bro0, bro1;
                cpline bound0, bound1;
                cpline edge0, edge1;
                cpoint* corner = nullptr, * pnt0 = nullptr, * pnt1 = nullptr;
                cpoint dir0, dir1;

                for (int i = 0; i < hexdata.size(); i++) {
                    memset(m_faults, 0, sizeof(int) * 8);

                    hexdata[i].cal_fault(&pdata[0], 1, m_faults);
                    get_ijk(i, cur_i, cur_j, cur_k, ti,tj,tk);

                    for (int fi = 0; fi < 8; fi++) {
                        if (m_faults[fi]!=0) {

                            corner = nullptr;
                            pnt0 = nullptr;
                            pnt1 = nullptr;
                            switch (fi)
                            {
                            case 0: {
                                bro0 = get_hex_id(ti, tj + 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti + 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[2]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[2]]);
                                corner = &(pdata[hexdata[i].mpoint[2]]);
                                dir0 = pdata[hexdata[bro0].mpoint[2]] - pdata[hexdata[bro0].mpoint[3]];
                                dir1 = pdata[hexdata[bro1].mpoint[2]] - pdata[hexdata[bro1].mpoint[1]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[2]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[1]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[3]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[1]], dir1, 1);
                                break;
                            }
                            case 1: {
                                bro0 = get_hex_id(ti, tj + 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti - 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[3]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[3]]);
                                corner = &(pdata[hexdata[i].mpoint[3]]);
                                dir0 = pdata[hexdata[bro0].mpoint[2]] - pdata[hexdata[bro0].mpoint[3]];
                                dir1 = pdata[hexdata[bro1].mpoint[0]] - pdata[hexdata[bro1].mpoint[3]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[2]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[0]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[2]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[0]], dir1, 1);
                                break;
                            }
                            case 2: {
                                bro0 = get_hex_id(ti, tj - 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti - 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[0]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[0]]);
                                corner = &(pdata[hexdata[i].mpoint[0]]);
                                dir0 = pdata[hexdata[bro0].mpoint[0]] - pdata[hexdata[bro0].mpoint[1]];
                                dir1 = pdata[hexdata[bro1].mpoint[0]] - pdata[hexdata[bro1].mpoint[3]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[0]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[0]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[1]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[3]], dir1, 1);
                                break;
                            }
                            case 3: {
                                bro0 = get_hex_id(ti, tj - 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti + 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[1]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[1]]);
                                corner = &(pdata[hexdata[i].mpoint[1]]);
                                dir0 = pdata[hexdata[bro0].mpoint[1]] - pdata[hexdata[bro0].mpoint[0]];
                                dir1 = pdata[hexdata[bro1].mpoint[1]] - pdata[hexdata[bro1].mpoint[2]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[1]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[1]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[0]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[2]], dir1, 1);
                                break;
                            }
                            case 4: {
                                bro0 = get_hex_id(ti, tj + 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti + 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[6]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[6]]);
                                corner = &(pdata[hexdata[i].mpoint[6]]);
                                dir0 = pdata[hexdata[bro0].mpoint[6]] - pdata[hexdata[bro0].mpoint[7]];
                                dir1 = pdata[hexdata[bro1].mpoint[6]] - pdata[hexdata[bro1].mpoint[5]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[6]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[6]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[7]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[5]], dir1, 1);
                                break;
                            }
                            case 5: {
                                bro0 = get_hex_id(ti, tj + 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti - 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[7]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[7]]);
                                corner = &(pdata[hexdata[i].mpoint[7]]);
                                dir0 = pdata[hexdata[bro0].mpoint[6]] - pdata[hexdata[bro0].mpoint[7]];
                                dir1 = pdata[hexdata[bro1].mpoint[4]] - pdata[hexdata[bro1].mpoint[7]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[7]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[7]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[6]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[4]], dir1, 1);
                                break;
                            }
                            case 6: {
                                bro0 = get_hex_id(ti, tj - 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti - 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0 == -1 || bro1 == -1 || !(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[4]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[4]]);
                                corner = &(pdata[hexdata[i].mpoint[4]]);
                                dir0 = pdata[hexdata[bro0].mpoint[4]] - pdata[hexdata[bro0].mpoint[5]];
                                dir1 = pdata[hexdata[bro1].mpoint[4]] - pdata[hexdata[bro1].mpoint[7]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[4]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[4]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[5]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[7]], dir1, 1);
                                break;
                            }
                            case 7: {
                                bro0 = get_hex_id(ti, tj - 1, tk, cur_i, cur_j, cur_k);
                                bro1 = get_hex_id(ti + 1, tj, tk, cur_i, cur_j, cur_k);
                                if (bro0==-1||bro1==-1||!(hexdata[bro0].stat == hexdata[bro1].stat && hexdata[bro0].stat == hexdata[i].stat))
                                    continue;
                                pnt0 = &(pdata[hexdata[bro0].mpoint[5]]);
                                pnt1 = &(pdata[hexdata[bro1].mpoint[5]]);
                                corner = &(pdata[hexdata[i].mpoint[5]]);
                                dir0 = pdata[hexdata[bro0].mpoint[5]] - pdata[hexdata[bro0].mpoint[4]];
                                dir1 = pdata[hexdata[bro1].mpoint[5]] - pdata[hexdata[bro1].mpoint[6]];
                                bound0 = cpline(pdata[hexdata[bro0].mpoint[5]], dir0, 1);
                                bound1 = cpline(pdata[hexdata[bro1].mpoint[5]], dir1, 1);
                                edge0 = cpline(pdata[hexdata[i].mpoint[4]], dir0, 1);
                                edge1 = cpline(pdata[hexdata[i].mpoint[6]], dir1, 1);
                                break;
                            }
                            }
                            if (corner != nullptr) {                            
                                *corner = edge0.intersection(edge1);
                                *pnt0 = bound0.intersection(edge1);
                                *pnt1 = bound1.intersection(edge0);
                            }
                        }
                    }
                }
            }
        }
        //更新储层ijk方向尺寸
        topo_i = cur_i;
        topo_j = cur_j;
        topo_k = cur_k;
        i_withext = topo_i;
        j_withext = topo_j;
        k_withext = topo_k;
        delete[]subhex;
    }

    //根据不同生长方向direct、生长层数layers、需求距离len来决定初始层高height、生长率h_ratio
    //direct表示生成方向，0是i方向，1是j方向，2是zup方向，3是zdown方向
    void confirm_height(int direct, double& height, double &h_ratio, int layers, int len) {
        double ref_h = 0.0;
        double* total_len;
        //获取cphex
        cphex t;
        t = hexdata[get_hex_id(
            reservior_relative_i_start,
            reservior_relative_j_start,
            reservior_relative_k_start,
            i_withext, j_withext, k_withext)];
        //匹配不同方向，获取参考高度和该方向的实际生长长度
        switch (direct)
        {
        case 0: {
            ref_h = fabs(pdata[t.mpoint[1]].x - pdata[t.mpoint[0]].x);
            total_len = &leni_ext;
            break;
        }
        case 1: {
            ref_h = fabs(pdata[t.mpoint[3]].y - pdata[t.mpoint[0]].y);
            total_len = &lenj_ext;

            break;
        }
        case 2:
            ref_h = fabs(pdata[t.mpoint[4]].z - pdata[t.mpoint[0]].z);
            total_len = &lenk_up_ext;
            break;
        case 3:
            ref_h = fabs(pdata[t.mpoint[4]].z - pdata[t.mpoint[0]].z);
            total_len = &lenk_bot_ext;
            break;
        default:
            break;
        }

        //将生长率h_ratio赋初值1.0
        double coef = 0, rat;
        h_ratio = 1.0;
        //迭代生长率，使之最终达到要求生长距离
        while (coef * ref_h < len) {
            coef = 0;
            rat = h_ratio;
            
            for (int i = 0; i <layers; i++) {
                coef += rat;
                rat *= h_ratio;
            }
            if (coef * ref_h >= len) {
                break;
            }
            
            if (coef*ref_h<len) {
                h_ratio += 0.01;
            }
            else {
                h_ratio -= 0.01;
            }
        }
        *total_len = coef * ref_h;

        height = h_ratio * ref_h;
    }

    //延拓网格
    //ratio_x/ratio_y:x、y方向网格延拓距离倍数
    //x_layers/y_layers/z_uplayers/z_bottomlayers:xyz方向网格延拓数量
    //z_positve_positon/z_negative_positon:z方向延拓距离
    void extent(double ratio_x, double ratio_y, int x_layers, int y_layers, int z_uplayers, int z_bottomlayers, double z_positve_positon, double z_negative_positon) {
        //内部数据存储的是海拔高度相反数
        swap(z_positve_positon, z_negative_positon);
        z_positve_positon *= -1.0;
        z_negative_positon *= -1.0;

        //根据模型包围盒，计算各方向需要延拓的长度
        double leni = (grid_bbox[1].x - grid_bbox[0].x) * ratio_x;
        double lenj = (grid_bbox[1].y - grid_bbox[0].y) * ratio_y;
        double lenz_up = z_positve_positon - grid_bbox[1].z;
        double lenz_bot = grid_bbox[0].z - z_negative_positon;

        //更新含延拓网格的ijk方向尺寸
        i_withext = topo_i + x_layers * 2;
        j_withext = topo_j + y_layers * 2;
        k_withext = topo_k + z_uplayers + z_bottomlayers;

        //更新储层ijk方向开始索引
        reservior_relative_i_start = x_layers;
        reservior_relative_j_start = y_layers;
        reservior_relative_k_start = z_bottomlayers;

        //更新储层ijk方向结束索引
        reservior_relative_i_end = reservior_relative_i_start + topo_i;
        reservior_relative_j_end = reservior_relative_j_start + topo_j;
        reservior_relative_k_end = reservior_relative_k_start + topo_k;

        //确认各方向首层高度和生长率
        confirm_height(0, i_ori_height, height_ratio[0], x_layers, leni);
        confirm_height(1, j_ori_height, height_ratio[1], y_layers, lenj);
        confirm_height(2, k_up_ori_height, height_ratio[2], z_uplayers, lenz_up);
        confirm_height(3, k_bot_ori_height, height_ratio[3], z_bottomlayers, lenz_bot);

        //将存储原始储层的hexdata扩容，同时并将原储层数据按原拓扑关系存储
        int hexid_withext = 0;
        vector<cphex> t_hexdata = hexdata;
        hexdata.clear();
        hexdata.resize(i_withext * j_withext * k_withext);
        int i_reservior = 0, j_reservior = 0, k_reservior = 0;
        for (int i = 0; i < t_hexdata.size(); i++) {
            get_ijk(i, topo_i, topo_j, topo_k, i_reservior, j_reservior, k_reservior);
            hexid_withext = get_hex_id(reservior_relative_i_start + i_reservior,
                reservior_relative_j_start + j_reservior,
                reservior_relative_k_start + k_reservior, i_withext, j_withext, k_withext);
            hexdata[hexid_withext] = t_hexdata[i];
        }

        //首先对六个面进行延拓
        cphex* host;
        cphex* extent_hex;
        int i_start, j_start, k_start, i_end, j_end, k_end;
        int i_offset = 0, j_offset = 0, k_offset = 0;
        int i_special_offset = 0;
        int j_special_offset = 0;
        int k_special_offset = 0;
        double rat = 1;

        //根据不同面，确定ijk方向起止索引，已经邻居的ijk坐标偏移值
        for (int j = reservior_relative_j_start; j < reservior_relative_j_end; j++) {
            for (int i = reservior_relative_i_start; i < reservior_relative_i_end; i++) {
                attach_hex_from_onehex(reservior_relative_k_start, 0,
                    i,
                    j,
                    reservior_relative_k_start);
                attach_hex_from_onehex(k_withext - reservior_relative_k_end, 5,
                    i,
                    j,
                    reservior_relative_k_end - 1);
            }
        }

        for (int k = reservior_relative_k_start; k < reservior_relative_k_end; k++) {
            for (int i = reservior_relative_i_start; i < reservior_relative_i_end; i++) {
                attach_hex_from_onehex(reservior_relative_j_start, 1,
                    i,
                    reservior_relative_j_start,
                    k);
                attach_hex_from_onehex(j_withext - reservior_relative_j_end, 3,
                    i,
                    reservior_relative_j_end - 1,
                    k);
            }
        }
        for (int j = reservior_relative_j_start; j < reservior_relative_j_end; j++) {
            for (int k = reservior_relative_k_start; k < reservior_relative_k_end; k++) {
                attach_hex_from_onehex(reservior_relative_i_start, 2,
                    reservior_relative_i_start,
                    j,
                    k);
                attach_hex_from_onehex(i_withext - reservior_relative_i_end, 4,
                    reservior_relative_i_end - 1,
                    j,
                    k);
            }
        }

        //沿上下底面的8条边方向进行延拓
        extent_edge_corner();
        //沿4条侧边进行延拓，同时也延拓8个角上位置的延拓网格
        //extent_corner();
        
        hexsize = hexdata.size();
        //标记延拓网格分区属性
        for (int k = 0; k < k_withext; ++k) {
            for (int i = 0; i < i_withext; i++) {
                for (int j = 0; j < j_withext; j++) {
                    if (k < reservior_relative_k_start) {
                        if (i >= reservior_relative_i_start && i < reservior_relative_i_end && j >= reservior_relative_j_start && j < reservior_relative_j_end) {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 4;
                        }
                        else {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 3;
                        }

                    }
                    else if (k < reservior_relative_k_end) {
                        if (i >= reservior_relative_i_start && i < reservior_relative_i_end && j >= reservior_relative_j_start && j < reservior_relative_j_end) {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 1;
                        }
                        else {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 2;
                        }
                    }
                    else {
                        if (i >= reservior_relative_i_start && i < reservior_relative_i_end && j >= reservior_relative_j_start && j < reservior_relative_j_end) {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 6;
                        }
                        else {
                            hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)].ext = 5;
                        }
                    }
                }
            }
        }
        return;
    }

    //延拓生长边角上的六面体单元
    void extent_edge_corner() {
        cphex* host1=nullptr;
        cphex* host2=nullptr;
        cphex* extent_hex=nullptr;
        //根据不同edge生长
        for (int j = reservior_relative_j_start; j < reservior_relative_j_end; ++j) {
            for (int i = reservior_relative_i_start-1; i >=0; i--) {
                for (int k = reservior_relative_k_start - 1; k >= 0; --k) {
                    attach_hex_from_2hex(i, j, k, 1);
                }
                for (int k = reservior_relative_k_end ; k <k_withext; ++k) {
                    attach_hex_from_2hex(i, j, k, 9);
                }
            }
        }
        for (int j = reservior_relative_j_start; j < reservior_relative_j_end; ++j) {
            for (int i = reservior_relative_i_end; i < i_withext; i++) {
                for (int k = reservior_relative_k_start - 1; k >= 0; --k) {
                    attach_hex_from_2hex(i, j, k, 3);
                }
                for (int k = reservior_relative_k_end; k < k_withext; ++k) {
                    attach_hex_from_2hex(i, j, k, 11);
                }
            }
        }

        for (int j = reservior_relative_j_start - 1; j >= 0; --j) {
            for (int i = reservior_relative_i_start; i < reservior_relative_i_end; i++) {
                for (int k = reservior_relative_k_start - 1; k >= 0; --k) {
                    attach_hex_from_2hex(i, j, k, 0);
                }
                for (int k = reservior_relative_k_end; k <k_withext; ++k) {
                    attach_hex_from_2hex(i, j, k, 8);
                }
            }
        }

        for (int j = reservior_relative_j_end; j < j_withext; ++j) {
            for (int i = reservior_relative_i_start; i < reservior_relative_i_end; i++) {
                for (int k = reservior_relative_k_start - 1; k >= 0; --k) {
                    attach_hex_from_2hex(i, j, k, 2);
                }
                for (int k = reservior_relative_k_end; k <k_withext; ++k) {
                    attach_hex_from_2hex(i, j, k, 10);
                }
            }
        }

        for (int k = 0; k < k_withext; ++k) {
            for (int j = reservior_relative_j_start - 1; j >= 0; --j) {
                for (int i = reservior_relative_i_start - 1; i >= 0; i--) {
                    attach_hex_from_2hex(i, j, k, 4);
                }
                for (int i = reservior_relative_i_end; i <i_withext; i++) {
                    attach_hex_from_2hex(i, j, k, 7);
                }
            }
            for (int j = reservior_relative_j_end; j < j_withext; ++j) {
                for (int i = reservior_relative_i_start - 1; i >= 0; i--) {
                    attach_hex_from_2hex(i, j, k, 5);
                }
                for (int i = reservior_relative_i_end; i < i_withext; i++) {
                    attach_hex_from_2hex(i, j, k, 6);
                }
            }
        }

        return;
    }

    void attach_hex_from_2hex(int i, int j, int k, int edge_id) {
        int edge_node[12][2] = {
            {0,1},{0,3},{3,2},{1,2},
            {0,4},{3,7},{2,6},{1,5},
            {4,5},{4,7},{7,6},{5,6}
        };
        //确定规则单元ijk长度
        double facet_len[3];
        facet_len[0] = (grid_bbox[1].x - grid_bbox[0].x) / topo_i;
        facet_len[1] = (grid_bbox[1].y - grid_bbox[0].y) / topo_j;
        facet_len[2] = (grid_bbox[1].z - grid_bbox[0].z) / topo_k;

        unordered_map<cpoint, int, cpoint_hash>::iterator it;
        cpoint surface_p1, surface_p2, host_p1, host_p2,target[2];
        cphex* host1, * host2, * sur, * temp;
        topo_hex h1;
        topo_hex h2;
        topo_hex ext;
        int ids[2];
        double rat = 0.0, len, surface;

        switch (edge_id)
        {
        case 0: {
            //根据不同方向确定生长方向上最外层的点
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k + 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j + 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i, j + 1, 0, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[0]].y;
            surface = pdata[sur->mpoint[0]].z;

            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(grid_bbox[0].x + (i - reservior_relative_i_start) * facet_len[0],
                len,
                surface);
            surface_p2 = surface_p1;
            surface_p2.x += facet_len[0];

            //将邻接两个单元与当前目标单元构建拓扑
            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXBOTTOM);
            ext.is_neighbor_of(h2, HEXFRONT);

            rat = fabs((pdata[temp->mpoint[7]].z - pdata[temp->mpoint[3]].z)) /
                fabs((pdata[sur->mpoint[2]].z - pdata[host1->mpoint[2]].z));

            //确定目标单元仍有的两个未知点
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            //点集去重
            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 1: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k + 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i + 1, j, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i + 1, j, 0, i_withext, j_withext, k_withext)];;
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[ids[0]]].x;
            surface = pdata[sur->mpoint[ids[0]]].z;

            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(len,
                grid_bbox[1].y - (j - reservior_relative_j_start) * facet_len[1],
                surface);
            surface_p2 = surface_p1;
            surface_p2.y -= facet_len[1];

            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXBOTTOM);
            ext.is_neighbor_of(h2, HEXLEFT);

            rat = fabs((pdata[temp->mpoint[6]].z - pdata[temp->mpoint[2]].z)) /
                fabs((pdata[sur->mpoint[ids[0]]].z - pdata[host1->mpoint[ids[0]]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 2: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k + 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j - 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i, j - 1, 0, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];
            len = pdata[host1->mpoint[3]].y;
            surface = pdata[sur->mpoint[0]].z;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(grid_bbox[0].x + (i - reservior_relative_i_start) * facet_len[0],
                len,
                surface);
            surface_p2 = surface_p1;
            surface_p2.x += facet_len[0];


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXBOTTOM);
            ext.is_neighbor_of(h2, HEXBACK);

            rat = fabs((pdata[temp->mpoint[4]].z - pdata[temp->mpoint[0]].z)) /
                fabs((pdata[sur->mpoint[0]].z - pdata[host1->mpoint[0]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 3: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k + 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i - 1, j, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i - 1, j, 0, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];
            len = pdata[host1->mpoint[1]].x;
            surface = pdata[sur->mpoint[0]].z;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(len,
                grid_bbox[1].y - (j - reservior_relative_j_start) * facet_len[1],
                surface);
            surface_p2 = surface_p1;
            surface_p2.y -= facet_len[1];


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXBOTTOM);
            ext.is_neighbor_of(h2, HEXRIGHT);

            rat = fabs((pdata[temp->mpoint[4]].z - pdata[temp->mpoint[0]].z)) /
                fabs((pdata[sur->mpoint[0]].z - pdata[host1->mpoint[0]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 4: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];

            host1 = &hexdata[get_hex_id(i + 1, j, k, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j + 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(0, j + 1, k, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[ids[0]]].z;
            surface = pdata[sur->mpoint[ids[0]]].x;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(pdata[sur->mpoint[ids[0]]].x,
                pdata[host1->mpoint[ids[0]]].y,
                pdata[sur->mpoint[ids[0]]].z);
            surface_p2 = cpoint(pdata[sur->mpoint[ids[1]]].x,
                pdata[host1->mpoint[ids[1]]].y,
                pdata[sur->mpoint[ids[1]]].z);


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXLEFT);
            ext.is_neighbor_of(h2, HEXFRONT);

            rat = fabs((pdata[temp->mpoint[2]].x - pdata[temp->mpoint[3]].x)) /
                fabs((pdata[sur->mpoint[0]].x - pdata[host1->mpoint[0]].x));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 5: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];

            host1 = &hexdata[get_hex_id(i + 1, j, k, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j - 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(0, j - 1, k, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(pdata[sur->mpoint[ids[0]]].x,
                pdata[host1->mpoint[ids[0]]].y,
                pdata[sur->mpoint[ids[0]]].z);
            surface_p2 = cpoint(pdata[sur->mpoint[ids[1]]].x,
                pdata[host1->mpoint[ids[1]]].y,
                pdata[sur->mpoint[ids[1]]].z);


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXLEFT);
            ext.is_neighbor_of(h2, HEXBACK);

            rat = fabs((pdata[temp->mpoint[0]].x - pdata[temp->mpoint[1]].x)) /
                fabs((pdata[sur->mpoint[0]].x - pdata[host1->mpoint[0]].x));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 6: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];

            host1 = &hexdata[get_hex_id(i - 1, j, k, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j - 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i_withext - 1, j - 1, k, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[ids[0]]].z;
            surface = pdata[sur->mpoint[ids[0]]].x;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(pdata[sur->mpoint[ids[0]]].x,
                pdata[host1->mpoint[ids[0]]].y,
                pdata[sur->mpoint[ids[0]]].z);
            surface_p2 = cpoint(pdata[sur->mpoint[ids[1]]].x,
                pdata[host1->mpoint[ids[1]]].y,
                pdata[sur->mpoint[ids[1]]].z);


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXRIGHT);
            ext.is_neighbor_of(h2, HEXBACK);

            rat = fabs((pdata[temp->mpoint[0]].x - pdata[temp->mpoint[1]].x)) /
                fabs((pdata[sur->mpoint[1]].x - pdata[host1->mpoint[1]].x));

            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 7: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];

            host1 = &hexdata[get_hex_id(i - 1, j, k, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j + 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i_withext - 1, j + 1, k, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[ids[0]]].z;
            surface = pdata[sur->mpoint[ids[0]]].x;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(pdata[sur->mpoint[ids[0]]].x,
                pdata[host1->mpoint[ids[0]]].y,
                pdata[sur->mpoint[ids[0]]].z);
            surface_p2 = cpoint(pdata[sur->mpoint[ids[1]]].x,
                pdata[host1->mpoint[ids[1]]].y,
                pdata[sur->mpoint[ids[1]]].z);


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXRIGHT);
            ext.is_neighbor_of(h2, HEXFRONT);

            rat = fabs((pdata[temp->mpoint[2]].x - pdata[temp->mpoint[3]].x)) /
                fabs((pdata[sur->mpoint[1]].x - pdata[host1->mpoint[1]].x));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 8: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k - 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j + 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i, j + 1, k_withext - 1, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];
            len = pdata[host1->mpoint[ids[0]]].y;
            surface = pdata[sur->mpoint[7]].z;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(grid_bbox[0].x + (i - reservior_relative_i_start) * facet_len[0],
                len,
                surface);
            surface_p2 = surface_p1;
            surface_p2.x += facet_len[0];


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXUP);
            ext.is_neighbor_of(h2, HEXFRONT);

            rat = fabs((pdata[temp->mpoint[7]].z - pdata[temp->mpoint[3]].z)) /
                fabs((pdata[sur->mpoint[7]].z - pdata[host1->mpoint[7]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 9: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k - 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i + 1, j, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i + 1, j, k_withext - 1, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];

            len = pdata[host1->mpoint[ids[0]]].x;
            surface = pdata[sur->mpoint[ids[0]]].z;

            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(len,
                grid_bbox[1].y - (j - reservior_relative_j_start) * facet_len[1],
                surface);
            surface_p2 = surface_p1;
            surface_p2.y -= facet_len[1];

            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXUP);
            ext.is_neighbor_of(h2, HEXLEFT);

            rat = fabs((pdata[temp->mpoint[6]].z - pdata[temp->mpoint[2]].z)) /
                fabs((pdata[sur->mpoint[ids[0]]].z - pdata[host1->mpoint[ids[0]]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 10: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k - 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i, j - 1, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i, j - 1, k_withext - 1, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];
            len = pdata[host1->mpoint[ids[0]]].y;
            surface = pdata[sur->mpoint[ids[0]]].z;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(grid_bbox[0].x + (i - reservior_relative_i_start) * facet_len[0],
                len,
                surface);
            surface_p2 = surface_p1;
            surface_p2.x += facet_len[0];


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXUP);
            ext.is_neighbor_of(h2, HEXBACK);

            rat = fabs((pdata[temp->mpoint[4]].z - pdata[temp->mpoint[0]].z)) /
                fabs((pdata[sur->mpoint[7]].z - pdata[host1->mpoint[7]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        case 11: {//见case0，代码功能基本一致
            ids[0] = edge_node[edge_id][0];
            ids[1] = edge_node[edge_id][1];
            host1 = &hexdata[get_hex_id(i, j, k - 1, i_withext, j_withext, k_withext)];
            host2 = &hexdata[get_hex_id(i - 1, j, k, i_withext, j_withext, k_withext)];

            sur = &hexdata[get_hex_id(i - 1, j, k_withext - 1, i_withext, j_withext, k_withext)];
            temp = &hexdata[get_hex_id(i, j, k, i_withext, j_withext, k_withext)];
            len = pdata[host1->mpoint[1]].x;
            surface = pdata[sur->mpoint[7]].z;
            host_p1 = pdata[host1->mpoint[ids[0]]];
            host_p2 = pdata[host1->mpoint[ids[1]]];

            surface_p1 = cpoint(len,
                grid_bbox[1].y - (j - reservior_relative_j_start) * facet_len[1],
                surface);
            surface_p2 = surface_p1;
            surface_p2.y -= facet_len[1];


            h1 = topo_hex(host1);
            h2 = topo_hex(host2);
            ext = topo_hex(temp);
            ext.is_neighbor_of(h1, HEXUP);
            ext.is_neighbor_of(h2, HEXRIGHT);

            rat = fabs((pdata[temp->mpoint[4]].z - pdata[temp->mpoint[0]].z)) /
                fabs((pdata[sur->mpoint[7]].z - pdata[host1->mpoint[7]].z));
            target[0] = (surface_p1 - host_p1) * rat + host_p1;
            target[1] = (surface_p2 - host_p2) * rat + host_p2;

            for (int ptr = 0; ptr < 2; ptr++) {
                it = point_map.find(target[ptr]);
                if (it == point_map.end()) {
                    point_map.insert({ target[ptr],pdata.size() });
                    pdata.push_back(target[ptr]);
                    *(ext.mpoint[ids[ptr]]) = pdata.size() - 1;;
                }
                else {
                    *(ext.mpoint[ids[ptr]]) = it->second;
                }
            }
            break;
        }
        default:
            break;
        }

    }


    // 0 bottom
    // 1 front
    // 2 left
    // 3 back
    // 4 right
    // 5 up
    //将延拓的最外层均当作是规则网格
    //将储层的最外层逐一与延拓最外层的规则网格对应相连，再根据首层高度、生长率、层数等划分出各层的节点，生成延拓网格
    void attach_hex_from_onehex(int layers,  int direct, int x, int y, int z) {
        double ratio = 1.0;
        cphex* host, * extent_hex;
        extent_hex = nullptr;
        host = &hexdata[get_hex_id(x, y, z, i_withext, j_withext, k_withext)];
        double facet_len[3];
        facet_len[0] = (grid_bbox[1].x - grid_bbox[0].x) / topo_i;
        facet_len[1] = (grid_bbox[1].y - grid_bbox[0].y) / topo_j;
        facet_len[2] = (grid_bbox[1].z - grid_bbox[0].z) / topo_k;

        //根据生长面方向以及ijk方向的索引x、y、z确定对应规则单元各点
        cpoint base_pnt;
        switch (direct) {
        case 0: {
            base_pnt = cpoint(grid_bbox[0].x + (x - reservior_relative_i_start) * facet_len[0],
                grid_bbox[1].y - (y - reservior_relative_j_start) * facet_len[1],
                grid_bbox[0].z - lenk_bot_ext);;
            break;
        }
        case 1: {
            base_pnt = cpoint(grid_bbox[0].x + (x - reservior_relative_i_start) * facet_len[0], 
                grid_bbox[1].y+lenj_ext, 
                grid_bbox[0].z + (z - reservior_relative_k_start) * facet_len[2]);
            break;
        }
        case 2: {
            base_pnt = cpoint(grid_bbox[0].x -leni_ext,
                grid_bbox[1].y - (y - reservior_relative_j_start) * facet_len[1],
                grid_bbox[0].z + (z - reservior_relative_k_start) * facet_len[2]);
            break;
        }
        case 3: {
            base_pnt = cpoint(grid_bbox[0].x + (x - reservior_relative_i_start) * facet_len[0],
                grid_bbox[0].y - lenj_ext + facet_len[1] ,
                grid_bbox[0].z + (z - reservior_relative_k_start) * facet_len[2]);
            break;
        }
        case 4: {
            base_pnt = cpoint(grid_bbox[1].x + leni_ext - facet_len[0], 
                grid_bbox[1].y - (y - reservior_relative_j_start) * facet_len[1],
                grid_bbox[0].z + (z - reservior_relative_k_start) * facet_len[2]);
            break;
        }
        case 5: {
            base_pnt = cpoint(grid_bbox[0].x + (x - reservior_relative_i_start) * facet_len[0],
                grid_bbox[1].y - (y - reservior_relative_j_start) * facet_len[1],
                grid_bbox[1].z + lenk_up_ext - facet_len[2]);
            break;
        }
        }
        cpoint base_hex[8] = { base_pnt, base_pnt, base_pnt, base_pnt, base_pnt, base_pnt, base_pnt, base_pnt };

        base_hex[1].x += facet_len[0];
        base_hex[2].x += facet_len[0];
        base_hex[2].y -= facet_len[1];
        base_hex[3].y -= facet_len[1];

        base_hex[4].z += facet_len[2];
        base_hex[5].x += facet_len[0];
        base_hex[5].z += facet_len[2];
        base_hex[6].x += facet_len[0];
        base_hex[6].y -= facet_len[1];
        base_hex[6].z += facet_len[2];
        base_hex[7].y -= facet_len[1];
        base_hex[7].z += facet_len[2];

        topo_hex to_host;
        topo_hex to_extent_hex;

        cpoint target[8];
        double rat = 1;
        double cal_len=0.0;
        double cur_rat = 0.0;
        cpoint ref_hex[8];
        cpoint t[4];
        
        //从储层最外层逐层向延拓网格最外层生长
        for (int i = 0; i < layers; i++) {
            if (direct == 0) {
                to_host = topo_hex(&hexdata[get_hex_id(x, y, z - i, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x , y, z - i - 1, i_withext, j_withext, k_withext)];
                cal_len += rat * k_bot_ori_height;
                rat *= height_ratio[3];
                cur_rat = cal_len / lenk_bot_ext;
            }
            else if (direct == 1) {
                to_host = topo_hex(&hexdata[get_hex_id(x, y - i, z, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x, y - i - 1, z, i_withext, j_withext, k_withext)];
                cal_len += rat * j_ori_height;
                rat *= height_ratio[1];
                cur_rat = cal_len / lenj_ext;
            }
            else if (direct == 2) {
                to_host = topo_hex(&hexdata[get_hex_id(x - i, y, z, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x - i - 1, y, z, i_withext, j_withext, k_withext)];
                cal_len += rat * i_ori_height;
                rat *= height_ratio[0];
                cur_rat = cal_len / leni_ext;
            }
            else if (direct == 3) {
                to_host = topo_hex(&hexdata[get_hex_id(x, y + i, z, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x, y + i + 1, z, i_withext, j_withext, k_withext)];
                cal_len += rat * j_ori_height;
                rat *= height_ratio[1];
                cur_rat = cal_len / lenj_ext;
            }
            else if (direct == 4) {
                to_host = topo_hex(&hexdata[get_hex_id(x + i, y, z, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x + i + 1, y, z, i_withext, j_withext, k_withext)];
                cal_len += rat * i_ori_height;
                rat *= height_ratio[0];
                cur_rat = cal_len / leni_ext;
            }
            else if (direct == 5) {
                to_host = topo_hex(&hexdata[get_hex_id(x, y, z + i, i_withext, j_withext, k_withext)]);
                extent_hex = &hexdata[get_hex_id(x, y, z + i + 1, i_withext, j_withext, k_withext)];
                cal_len += rat * k_up_ori_height;
                rat *= height_ratio[2];
                cur_rat = cal_len / lenk_up_ext;
            }
            
            to_extent_hex = topo_hex(extent_hex);
            to_extent_hex.is_neighbor_of(to_host, direct);

            for (int m = 0; m < 8; m++) {
                target[m] = (base_hex[m] - pdata[host->mpoint[m]]) * cur_rat + pdata[host->mpoint[m]];
            }

            int face_pnt_id[6][4]{
                {0,1,2,3},
                {0,1,4,5},
                {0,3,4,7},
                {3,2,7,6},
                {1,2,5,6},
                {4,5,6,7} };

            int pid = 0;
            auto it = point_map.find(target[0]);
            for (int i = 0; i < 4; i++) {
                pid = face_pnt_id[direct][i];
                it = point_map.find(target[pid]);
                if (it == point_map.end()) {
                    point_map.insert({ target[pid], pdata.size() });
                    pdata.push_back(target[pid]);
                    extent_hex->mpoint[pid] = pdata.size() - 1;
                }
                else {
                    extent_hex->mpoint[pid] = it->second;
                }
            }
        }
    }

    //获取原始网格的子单元id，输入参数为在CPgrid中的ijk尺寸
    vector<int> children_id_of(int root_i,int root_j,int root_k){
        if (root_i >= 0 && root_i < ori_ijk[0] &&
            root_j >= 0 && root_j < ori_ijk[1] &&
            root_k >= 0 && root_k < ori_ijk[2]) {
            int cphex_times = 1;
            for (int i = 0; i < ref_times; i++) {
                cphex_times *= 2;
            }
            vector<int> ids;
            int istart = root_i * cphex_times;
            int jstart = root_j * cphex_times;
            int kstart = root_k * cphex_times;

            int iend = (root_i + 1) * cphex_times;
            int jend = (root_j + 1) * cphex_times;
            int kend = (root_k + 1) * cphex_times;

            for (int i = istart; i < iend; i++) {
                for (int j = jstart; j < jend; j++) {
                    for (int k = kstart; k < kend; k++) {
                        ids.push_back(get_hex_id(i + reservior_relative_i_start, j + reservior_relative_j_start, k + reservior_relative_k_start, i_withext, j_withext, k_withext));
                    }
                }
            }
            return ids;
        }
        else {
            cout << "Wrong index of root!:children_id_of\n";
            return vector<int> {};
        }
    }

    //获取finiteGrid中子单元的原始单元
    //输入参数为finiteGrid中子单元的ijk
    vector<int> root_ijk_of(int i,int j,int k){
        if (i < reservior_relative_i_start || i >= reservior_relative_i_end ||
            j < reservior_relative_j_start || j >= reservior_relative_j_end ||
            k < reservior_relative_k_start || k >= reservior_relative_k_end) {
            cout << "index out of reservior in:root_ijk_of\n";
            return vector<int>{};
        }

        int cphex_times=1;
        for(int i=0;i<ref_times;i++){
            cphex_times*=2;
        }

        i -= reservior_relative_i_start;
        j -= reservior_relative_j_start;
        k -= reservior_relative_k_start;
        
        i/=cphex_times;
        j/=cphex_times;
        k/=cphex_times;

        return vector<int>{i, j, k};
    }

    int add_cpoint(cphex* extent_hex, cpoint target[8],int direct) {
        // 0 bottom
        // 1 front
        // 2 left
        // 3 back
        // 4 right
        // 5 up
        int face_pnt_id[6][4]{
            {0,1,2,3},
            {0,1,4,5},
            {0,3,4,7},
            {3,2,7,8},
            {1,2,5,6},
            {4,5,6,7} };

        int pid = 0;
        auto it = point_map.find(target[0]);
        for (int i = 0; i < 4; i++) {
            pid = face_pnt_id[direct][i];
            it= point_map.find(target[pid]);

            if (it == point_map.end()) {
                point_map.insert({ target[pid], pdata.size() });
                pdata.push_back(target[pid]);
                extent_hex->mpoint[i] = pdata.size() - 1;
            }
            else {
                extent_hex->mpoint[i] = it->second;
            }
        }
    }

    void get_ijk(int id, int i, int j, int k, int& x, int& y, int& z) {
        x = id % i;
        y = id % (i * j) / i;
        z = id / (i * j);
    }

    int get_hex_id(int x, int y, int z, int i, int j, int k) {
        if (x >= 0 && x < i &&
            y >= 0 && y < j &&
            z >= 0 && z < k)
        {
            return x + y * i + z * i * j;
        }
        else
        {
            return -1;
        }
    }

    void optimize() {
        for (int i = 0; i < hexdata.size(); i++) {
            hexdata[i].optimize(pdata, point_map);
        }
    }

};
#endif